
clear all
set more off
capture log close

global dta "D:\Work\Valuing public good WTP\dta"
global results "D:\Work\Valuing public good WTP\Results"

*500 iterations - Oct 2021
*****************************
*  Probability as function  *
*****************************

foreach n in 1 2 3 4 5 6 7 8 9 10{
gen varname`n'=.
}
drop if _n>=1
save "$dta/MonteCarlo_Results_joint.dta",replace

  program define dh_dc
    version 14
    args lnf beta1 beta2 beta6 beta3 beta4 beta5 beta7
    tempvar d non_l pz p0 p1
    quietly gen double `d'=($ML_y1==1)

	quietly gen double `non_l'=(`beta3'-normal(`beta4'+`beta5')*`beta6')/(normal(`beta4'+`beta5')-normal(`beta4'))

	quietly gen double `pz'=binormal(`beta1',`beta2'+`non_l',`beta7')
    quietly gen double `p0'=1-(`pz')
    quietly gen double `p1'=(`pz')
    quietly replace `lnf'=(1-`d')*ln(`p0')+`d'*ln(`p1')
 end 

foreach N in 1000 3000 5000 10000 20000 { 
clear
set seed 123456789
forval i = 1/500{
clear
set obs `N'
di "i " `i'
**********
* Setup  *
**********
*Generate error terms
gen ee = -log(-log(runiform()))

*generate joint standard normal error terms
scalar rho=.5
matrix C=(1,rho \ rho,1)
drawnorm ee1 e, cov(C)
sum ee1
sum e
reg ee1 e

*Generate participation variables
gen Age=runiform(18,70)
gen Inc=exp(rnormal(4,1))

*Generate participation parameters
scalar gama_Age=0.003
scalar gama_Inc=.0014
scalar gama_C=-1

*Participation decision and prob 
gen Participate=gama_C+gama_Age*Age+gama_Inc*Inc+e

*Generate donation variables
gen A=exp(rnormal(0,2))

*Generate utility parameters
scalar alfa_A=.4
scalar alfa_Age=.05
scalar alfa_Inc=.027
scalar beta=.1

*------- Senario 1-model subjective probability as a function
*Generate subjective probability variables

gen D=runiform(1000,5000)
scalar D_threshold=D
drop D
di D_threshold
gen D_bar=D_threshold
gen c_solicAmount=floor(runiform(5,120))
*Generate probability parameters
scalar miu=0
scalar tao=.02
scalar ro=-.1/`N'
*Generate probability as functions
*
gen P_D1=normprob(miu+tao*c_solicAmount+ro*D_threshold)
gen P_D0=normprob(miu+ro*D_threshold)

*Now simulate the decisions 
*Generate participation decision
gen Part=(Participate>0)
gen WG=1
scalar wg=WG
gen Delta_EU = WG+(P_D1-P_D0)*(alfa_A*A+alfa_Age*Age+alfa_Inc*Inc+ee1)-beta*P_D1*c_solicAmount
gen Di=(Part==1&Delta_EU>0)
tab Part
tab Di

************
* Estimate *
************
 global varlist_ydc "A Age Inc "
 global varlist_d "Age Inc "

 global price "c_solicAmount"
 sum $varlist_d

 probit Part $varlist_d
 mat list e(b)
 mat bparticipation=e(b)
 mat list bparticipation

 scalar mu_prime=miu+ro*D_threshold
 mat b=(bparticipation,alfa_A,alfa_Age,alfa_Inc,beta,wg,mu_prime,tao,rho)

 ml model lf dh_dc (hurdle: Di=$varlist_d)(above: Di=$varlist_ydc,nocons)(beta:$price,nocons)(warm_glow:)(mu:)(tao:$price,nocons)(rho:)
 ml init b,copy
 
 cap {
 ml maximize,iterate(500) dif
 estat ic
 
 mat list e(b)
 di _b[above: A]
 di _b[above: Age]
 di _b[above: Inc]
 di _b[rho:_cons]
 }
 cap {
 mat V=e(V)
 di V[1,1]
 gen SEsq_A=V[4,4]
 gen SEsq_Age=V[5,5]
 gen SEsq_Inc=V[6,6]
 gen est_A=(V[4,4]!=0)
 gen est_Age=(V[5,5]!=0)
 gen est_Inc=(V[6,6]!=0)
 
 gen WTP_est=(_b[above: A]*A*est_A+_b[above: Age]*Age*est_Age+_b[above: Inc]*Inc*est_Inc)/_b[beta:c_solicAmount]
 sum WTP_est

 scalar WTP_e_scalar=r(mean)
 scalar WTP_sd_scalar=r(sd)
 scalar WTP_95Ci_lb=r(mean)-r(sd)*1.96
 scalar WTP_95Ci_ub=r(mean)+r(sd)*1.96
 scalar WTP_90Ci_lb=r(mean)-r(sd)*1.645
 scalar WTP_90Ci_ub=r(mean)+r(sd)*1.645
 
  *Generate the reference value of WTP here
 gen WTP=(.4*A+.05*Age+.027*Inc)/.1
 sum WTP
 scalar WTP_o=r(mean)
 scalar delta_WTP = WTP_e_scalar-WTP_o
 scalar total_benefit=e(N) * WTP_o
 scalar benefit_cost_ratio=total_benefit/D_threshold
 
 di WTP_95Ci_lb
 di WTP_95Ci_ub
 gen CI95_cover=1 if WTP_95Ci_ub>WTP_o & WTP_o>WTP_95Ci_lb
 scalar CI95cover=CI95_cover
 gen CI90_cover=1 if WTP_90Ci_ub>WTP_o & WTP_o>WTP_90Ci_lb
 scalar CI90cover=CI90_cover
 count if Di==1
 scalar N_donor=r(N)
 
 mat R=(`i',delta_WTP, WTP_e_scalar, WTP_o, CI90cover, CI95cover, e(N), total_benefit, D_threshold, benefit_cost_ratio,_b[beta:c_solicAmount],_b[above: A],_b[above: Age],_b[above: Inc],_b[tao:c_solicAmount],N_donor)
 mat list R
 
 clear 
 tempfile R_temp
 svmat R, names(varname)
 save "$dta/temp_result.dta",replace
 
 use "$dta/MonteCarlo_Results_joint.dta",clear
 append using "$dta/temp_result.dta"
 save "$dta/MonteCarlo_Results_joint.dta",replace
 }
 }
 }

 
 

************Compute Error Measurements***************
*Scenario 1-model subjective probability as a variable
foreach N in 1000 3000 5000 10000 20000  { 
clear all
use "$dta/MonteCarlo_Results_joint.dta",clear
quietly {
ren varname1 ID
ren varname2 delta_WTP
ren varname3 WTP_est
ren varname4 WTP_obs
ren varname5 CI90cover
ren varname6 CI95cover
ren varname7 N
ren varname8 total_benefit
ren varname9 total_cost
ren varname10 benefit_cost_ratio
ren varname11 beta
ren varname12 alfa_A
ren varname13 alfa_age
ren varname14 alfa_inc
ren varname15 tao_inc_probpara
ren varname16 Ndonors

keep if N==`N'
count if delta_WTP>1000|delta_WTP<-1000
*12 outliners not count
drop if delta_WTP>1000|delta_WTP<-1000

gen order=_n
egen success_N=max(order)
egen total_N=max(ID)
scalar Success_rate=success_N/total_N
di Success_rate

replace delta_WTP=WTP_est-WTP_obs
sum delta_WTP
scalar Mean_bias=r(mean)
di Mean_bias

sum delta_WTP,d
scalar Median_bias=r(p50)
di Median_bias 

gen sq_delta_WTP=delta_WTP^2
sum sq_delta_WTP 
scalar RMSE=sqrt(r(mean))
di RMSE

gen abs_delta_WTP=abs(delta_WTP)
sum abs_delta_WTP,d
scalar MAE=sqrt(r(p50))
di MAE

sum WTP_est
scalar MEAN_est=r(mean)
scalar SD_est=r(sd)
di MEAN_est " " SD_est

sum WTP_obs
scalar MEAN_WTP=r(mean)
di MEAN_WTP

egen CI90_covers=sum(CI90cover)
scalar CI90_coverrate=CI90_covers/success_N
di CI90_coverrate

egen CI95_covers=sum(CI95cover)
scalar CI95_coverrate=CI95_covers/success_N
di CI95_coverrate

sum Ndonors
scalar MEAN_donors=r(mean)
scalar SD_donors=r(sd)
di MEAN_donors
di SD_donors

sum beta
scalar MEAN_beta=r(mean)
scalar SD_beta=r(sd)

sum alfa_A
scalar MEAN_alfa_A=r(mean)
scalar SD_alfa_A=r(sd)

sum alfa_age
scalar MEAN_alfa_age=r(mean)
scalar SD_alfa_age=r(sd)

sum alfa_inc
scalar MEAN_alfa_inc=r(mean)
scalar SD_alfa_inc=r(sd)
}
di "N=" `N'
di "Success_rate " Success_rate " MEAN_est " MEAN_est " SD_est " SD_est " MEAN_WTP " MEAN_WTP   

di "Mean_bias " Mean_bias " Median_bias " Median_bias " RMSE " RMSE " MAE " MAE " CI90_coverrate " CI90_coverrate " CI95_coverrate " CI95_coverrate

eststo clear
estadd scalar Samplesize=`N'
estadd scalar Success_rate=Success_rate
estadd scalar MEAN_est=MEAN_est
estadd scalar SD_est=SD_est
estadd scalar MEAN_WTP=MEAN_WTP
estadd scalar Mean_bias=Mean_bias
estadd scalar Median_bias=Median_bias
estadd scalar RMSE=RMSE
estadd scalar MAE=MAE
estadd scalar CI90_coverrate=CI90_coverrate
estadd scalar CI95_coverrate=CI95_coverrate
estadd scalar MEAN_donors=MEAN_donors
estadd scalar SD_donors=SD_donors
estadd scalar MEAN_beta=MEAN_beta
estadd scalar SD_beta=SD_beta
estadd scalar MEAN_alfa_A=MEAN_alfa_A
estadd scalar SD_alfa_A=SD_alfa_A
estadd scalar MEAN_alfa_age=MEAN_alfa_age
estadd scalar SD_alfa_age=SD_alfa_age
estadd scalar MEAN_alfa_inc=MEAN_alfa_inc
estadd scalar SD_alfa_inc=SD_alfa_inc

esttab using"$results\Simu_results_joint_`N'.csv", replace b(a3) mti("") stats (Samplesize Success_rate MEAN_est SD_est MEAN_WTP Mean_bias Median_bias RMSE MAE CI90_coverrate CI95_coverrate MEAN_donors SD_donors MEAN_beta SD_beta MEAN_alfa_A SD_alfa_A MEAN_alfa_age SD_alfa_age MEAN_alfa_inc SD_alfa_inc)
}






*******************************
*  Probability as a variable  *
*******************************
*--------Scenario 2-model subjective probability as a variable
*It's the same P_D1 and P_D0, just that they are known to researchers in the estimation
clear all
set more off
global dta "D:\Work\Valuing public good WTP\dta"

foreach n in 1 2 3 4 5 6 7 8 9 10{
gen varname`n'=.
}
drop if _n>=1
save "$dta/MonteCarlo_Results1.dta",replace

global P_D1 "P_D1"
global P_D0 "P_D0"

   program define dh_dc_2
    version 14
    args lnf beta1 beta2 beta6 beta3 beta7
    tempvar d non_l pz p0 p1
    quietly gen double `d'=($ML_y1==1)

	quietly gen double `non_l'=(`beta3'-$P_D1*`beta6')/($P_D1-$P_D0)

	quietly gen double `pz'=binormal(`beta1',`beta2'+`non_l',`beta7')
    quietly gen double `p0'=1-(`pz')
    quietly gen double `p1'=(`pz')
    quietly replace `lnf'=(1-`d')*ln(`p0')+`d'*ln(`p1')
 end 
 
foreach N in 1000 3000 5000 10000 20000 { 
clear
set seed 12345678

forval i = 1/500{
clear
set obs `N'
di "i " `i'
**********
* Setup  *
**********
*Generate error terms
drawnorm e
gen ee = -log(-log(runiform()))
drawnorm ee1
gen ee2 = invnorm(runiform())

*Generate participation variables
gen Age=runiform(18,70)
gen Inc=exp(rnormal(4,1))
*Generate participation parameters
scalar gama_Age=.003
scalar gama_Inc=.0014
scalar gama_C=-1

*Participation decision and prob 
gen Participate=gama_C+gama_Age*Age+gama_Inc*Inc+e

*Generate donation variables
gen A=exp(rnormal(0,2))
*Generate utility parameters
scalar alfa_A=.4
scalar alfa_Age=.05
scalar alfa_Inc=.027
scalar beta=.1

*------- Senario 2
*Generate subjective probability variables
gen epsilon=rnormal(0,.1)
gen D=runiform(1000,5000)
scalar D_threshold=D
drop D
di D_threshold
gen D_bar=D_threshold
gen c_solicAmount=floor(runiform(5,120))
*Generate probability parameters
scalar miu=0
scalar tao=.02
scalar ro=-.1/`N'
*Generate probability as functions
gen P_D1=normprob(miu+tao*c_solicAmount+ro*D_threshold)
gen P_D0=normprob(miu+ro*D_threshold)
*Notice: though they are generated from functions, the probabilities will serve as observed variables in the estimation.

*Now simulate the decisions 
*Generate participation decision
gen Part=(Participate>0)
gen WG=1
scalar wg=WG
*Generate expected utility comparison
gen Delta_EU = 1+(P_D1-P_D0)*(alfa_A*A+alfa_Age*Age+alfa_Inc*Inc+ee1)-beta*P_D1*c_solicAmount
gen Di=(Part==1&Delta_EU>0)
tab Di

************
* Estimate *
************
 *  B C Education
 global varlist_ydc "A Age Inc "
 * Education
 global varlist_d "Age Inc "

 global price "c_solicAmount"
 sum $varlist_d

 probit Part $varlist_d
 mat list e(b)
 mat bparticipation=e(b)

 mat b=(bparticipation,alfa_A,alfa_Age,alfa_Inc,beta,wg)
 *D_bar
 ml model lf dh_dc_2 (hurdle: Di=$varlist_d)(above: Di=$varlist_ydc,nocons)(beta:$price,nocons)(warm_glow:)(rho:)
 ml init b,copy
 
 cap {
 ml maximize,iterate(500) dif
 estat ic

 mat list e(b)
 di _b[above: A]
 di _b[above: Age]
 di _b[above: Inc]
 di _b[rho:_cons]
 }
  
 cap{
 mat V=e(V)
 di V[1,1]
 gen SEsq_A=V[4,4]
 gen SEsq_Age=V[5,5]
 gen SEsq_Inc=V[6,6]
 gen est_A=(V[4,4]!=0)
 gen est_Age=(V[5,5]!=0)
 gen est_Inc=(V[6,6]!=0)
 
 gen WTP_est=(_b[above: A]*A*est_A+_b[above: Age]*Age*est_Age+_b[above: Inc]*Inc*est_Inc)/_b[beta:c_solicAmount]
 sum WTP_est
 scalar WTP_e_scalar=r(mean)
 scalar WTP_sd_scalar=r(sd)
 scalar WTP_95Ci_lb=r(mean)-r(sd)*1.96
 scalar WTP_95Ci_ub=r(mean)+r(sd)*1.96
 scalar WTP_90Ci_lb=r(mean)-r(sd)*1.645
 scalar WTP_90Ci_ub=r(mean)+r(sd)*1.645
 
 gen WTP=(.4*A+.05*Age+.027*Inc)/.1
 sum WTP
 scalar WTP_o=r(mean)
 scalar delta_WTP = WTP_e_scalar-WTP_o
 scalar total_benefit=e(N) * WTP_o
 scalar benefit_cost_ratio=total_benefit/D_threshold
 
 di WTP_95Ci_lb
 di WTP_95Ci_ub
 gen CI95_cover=1 if WTP_95Ci_ub>WTP_o & WTP_o>WTP_95Ci_lb
 scalar CI95cover=CI95_cover
 gen CI90_cover=1 if WTP_90Ci_ub>WTP_o & WTP_o>WTP_90Ci_lb
 scalar CI90cover=CI90_cover
 count if Di==1
 scalar N_donor=r(N)
 
 mat R=(`i',delta_WTP, WTP_e_scalar, WTP_o, CI90cover, CI95cover, e(N), total_benefit, D_threshold, benefit_cost_ratio,_b[beta:c_solicAmount],_b[above: A],_b[above: Age],_b[above: Inc],N_donor)
 mat list R
 
 clear 
 tempfile R_temp
 svmat R, names(varname)
 save "$dta/temp_result.dta",replace
 
 use "$dta/MonteCarlo_Results1.dta",clear
 append using "$dta/temp_result.dta"
 save "$dta/MonteCarlo_Results1.dta",replace
 }
 }
 }
*Scenario 2-model subjective probability as a variable
*50000
foreach N in 1000 3000 5000 10000 20000 { 
clear all
use "$dta/MonteCarlo_Results1.dta",clear
quietly {
ren varname1 ID
ren varname2 delta_WTP
ren varname3 WTP_est
ren varname4 WTP_obs
ren varname5 CI90cover
ren varname6 CI95cover
ren varname7 N
ren varname8 total_benefit
ren varname9 total_cost
ren varname10 benefit_cost_ratio
ren varname11 beta
ren varname12 alfa_A
ren varname13 alfa_age
ren varname14 alfa_inc
ren varname15 Ndonors

keep if N==`N'
count if delta_WTP>1000|delta_WTP<-1000
*12 outliners not count
drop if delta_WTP>1000|delta_WTP<-1000

gen order=_n
egen success_N=max(order)
egen total_N=max(ID)
scalar Success_rate=success_N/total_N
di Success_rate

replace delta_WTP=WTP_est-WTP_obs
sum delta_WTP
scalar Mean_bias=r(mean)
di Mean_bias

sum delta_WTP,d
scalar Median_bias=r(p50)
di Median_bias 

gen sq_delta_WTP=delta_WTP^2
sum sq_delta_WTP 
scalar RMSE=sqrt(r(mean))
di RMSE

gen abs_delta_WTP=abs(delta_WTP)
sum abs_delta_WTP,d
scalar MAE=sqrt(r(p50))
di MAE

sum WTP_est
scalar MEAN_est=r(mean)
scalar SD_est=r(sd)
di MEAN_est " " SD_est

sum WTP_obs
scalar MEAN_WTP=r(mean)
di MEAN_WTP

egen CI90_covers=sum(CI90cover)
scalar CI90_coverrate=CI90_covers/success_N
di CI90_coverrate

egen CI95_covers=sum(CI95cover)
scalar CI95_coverrate=CI95_covers/success_N
di CI95_coverrate

sum Ndonors
scalar MEAN_donors=r(mean)
scalar SD_donors=r(sd)
di MEAN_donors
di SD_donors

sum beta
scalar MEAN_beta=r(mean)
scalar SD_beta=r(sd)

sum alfa_A
scalar MEAN_alfa_A=r(mean)
scalar SD_alfa_A=r(sd)

sum alfa_age
scalar MEAN_alfa_age=r(mean)
scalar SD_alfa_age=r(sd)

sum alfa_inc
scalar MEAN_alfa_inc=r(mean)
scalar SD_alfa_inc=r(sd)
}
di "N=" `N'
di "Success_rate " Success_rate " MEAN_est " MEAN_est " SD_est " SD_est " MEAN_WTP " MEAN_WTP   

di "Mean_bias " Mean_bias " Median_bias " Median_bias " RMSE " RMSE " MAE " MAE " CI90_coverrate " CI90_coverrate " CI95_coverrate " CI95_coverrate
eststo clear
estadd scalar Samplesize=`N'
estadd scalar Success_rate=Success_rate
estadd scalar MEAN_est=MEAN_est
estadd scalar SD_est=SD_est
estadd scalar MEAN_WTP=MEAN_WTP
estadd scalar Mean_bias=Mean_bias
estadd scalar Median_bias=Median_bias
estadd scalar RMSE=RMSE
estadd scalar MAE=MAE
estadd scalar CI90_coverrate=CI90_coverrate
estadd scalar CI95_coverrate=CI95_coverrate
estadd scalar MEAN_donors=MEAN_donors
estadd scalar SD_donors=SD_donors
estadd scalar MEAN_beta=MEAN_beta
estadd scalar SD_beta=SD_beta
estadd scalar MEAN_alfa_A=MEAN_alfa_A
estadd scalar SD_alfa_A=SD_alfa_A
estadd scalar MEAN_alfa_age=MEAN_alfa_age
estadd scalar SD_alfa_age=SD_alfa_age
estadd scalar MEAN_alfa_inc=MEAN_alfa_inc
estadd scalar SD_alfa_inc=SD_alfa_inc

esttab using"$results\Simu_results1_5000inter_`N'.csv", replace b(a3) mti("") stats (Samplesize Success_rate MEAN_est SD_est MEAN_WTP Mean_bias Median_bias RMSE MAE CI90_coverrate CI95_coverrate MEAN_donors SD_donors MEAN_beta SD_beta MEAN_alfa_A SD_alfa_A MEAN_alfa_age SD_alfa_age MEAN_alfa_inc SD_alfa_inc)

}
